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SPATIAL MODELING OF EXTREME SNOW DEPTH 

By Juliette Blanchet and Anthony C. Davison^ 

Ecole Polytechnique Federate de Lausanne 

The spatial modeling of extreme snow is important for adequate 
risk management in Alpine and high altitude countries. A natu- 
ral approach to such modeling is through the theory of max-stable 
processes, an infinite-dimensional extension of multivariate extreme 
value theory. In this paper we describe the application of such pro- 
cesses in modeling the spatial dependence of extreme snow depth 
in Switzerland, based on data for the winters 1966-2008 at 101 sta- 
tions. The models we propose rely on a climate transformation that 
allows us to account for the presence of climate regions and for direc- 
tional effects, resulting from synoptic weather patterns. Estimation is 
performed through pairwise likelihood inference and the models are 
compared using penalized likelihood criteria. The max-stable models 
provide a much better fit to the joint behavior of the extremes than 
do independence or full dependence models. 

1. Introduction. Heavy snow events are among the most severe natu- 
ral hazards in mountainous countries. Every year, winter storms can hinder 
mobility by disrupting rail, road and air traffic. Extreme snowfall can over- 
load buildings and cause them to collapse, and can lead to flooding due 
to subsequent melting. Deep snow, combined with strong winds and un- 
stable snowpack, contributes to the formation of avalanches, and can cause 
fatalities and economic loss due to property damage or reduced mobility. 
The quantitative analysis of extreme snow events is important for the di- 
mensioning of avalanche defence structures, bridges and buildings, for flood 
protection measures and for integral risk management. 

Compared to phenomena such as rain, wind or temperature, extreme- 
value statistics of snow has been little studied. Bocchiola, Medagliani and 
Rosso (2006) and Bocchiola et al. (2008) analyzed three-day snowfall depth 
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in the Italian and Swiss Alps, and more recently Blancliet, Marty and Lelin- 
ing (2009) analyzed extreme snowfall in Switzerland. These articles derive 
characteristics of extreme snow events based on univariate extreme-value 
modeling which does not account for the dependence across different sta- 
tions. The spatial dependence of extreme snow data has yet to be discussed 
in the literature. 

Statistical modeling with multivariate extreme value distributions be- 
gan around two decades ago with publications such as Tawn (1988) and 
Coles and Tawn (1991), and has subsequently often been used for quan- 
tifying extremal dependence in applications. Financial examples are cur- 
rency exchange rate data [Hauksson et al. (2001)], swap rate data [Hsing, 
Kliippelberg and Kuhn (2004)] and stock market returns [Poon, Rockinger 
and Tawn (2003, 2004)], and environmental examples are rainfall data [Schlat- 
her and Tawn (2003)], oceanographical data [de Haan and de Ronde (1998); 
Coles and Tawn (1994)] and wind speed data [Coles and Walshaw (1994); 
Fawcett and Walshaw (2006)]. None of these articles treats the process under 
study as a spatial extension of multivariate extreme value theory. 

Until recently, a key difficulty in studying extreme events of spatial pro- 
cesses has been the lack of flexible models and appropriate inferential tools. 
Two different approaches to overcome this have been proposed. The first and 
most popular is to introduce a latent process, conditional on which standard 
extreme models are applied [Coles and Casson (1998); Fawcett and Walshaw 
(2006); Cooley et al. (2006); Cooley, Nychka and Naveau (2007); Gaetan and 
Grigoletto (2007); Sang and Gelfand (2009b); Eastoe (2009)]. Such models 
can be fitted using Markov chain Monte Carlo simulation, but they postu- 
late independence of extremes conditional on the latent process, and this 
is implausible in applications. One approach to introducing dependence is 
through a spatial copula, as suggested by Sang and Gelfand (2009a), but 
although this approach is an improvement, Davison, Padoan and Ribatet 
(2010) show that it can nevertheless lead to inadequate modeling of ex- 
treme rainfall. A second approach now receiving increasing attention rests 
on max-stable processes, first suggested by de Haan (1984) and developed 
by, for example, Schlather (2002) and Kabluchko, Schlather and de Haan 
(2009). Recent applications to rainfall data can be found in Buishand, de 
Haan and Zhou (2008), Smith and Stephenson (2009), Padoan, Ribatet and 
Sisson (2010) and Davison, Padoan and Ribatet (2010), and to temperature 
data in Davison and Gholamrezaee (2010). Max-stable modeling has the po- 
tential advantage of accounting for spatial dependence of extremes in a way 
that is consistent with the classical extreme-value theory, but is much less 
well developed than the use of latent processes or copulas. 

In the present paper, we use data from a denser measurement network 
than for previous applications. Owing to complex topography and weather 
patterns, the processes of Schlather (2002) and Smith (1990) cannot account 
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for the joint distribution of the extremes, and we therefore propose more 
complex models. We begin with an exploratory analysis highlighting some 
of the peculiarities of the data, and then in Section 3 present the max- 
stable processes of Schlather (2002) and Smith (1990), which are extended 
in Section 4 to our extreme snow depth data. As full likelihood inference is 
impossible for such models, in Section 5 we discuss how composite likelihood 
inference may be used for model estimation and comparison. The results of 
the data analysis are presented in Section 6 and a concluding discussion is 
given in Section 7. 

2. Preliminaries. 



2.1. Data. We consider annual maximum snow depth from the 101 sta- 
tions whose locations are shown in Figure 1. The stations belong to two 
networks run by the WSL Institute for Snow and Avalanche Research (SLF) 
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Fig. 1. Topography and locations of stations for which daily snow depth data are avail- 
able. First row: Topographical map of Switzerland (left) and station locations (right). Sec- 
ond row: Histogram of elevation of Switzerland at a 1 km grid (left) and of the stations 
(right). Color indicates altitude in meters above mean sea level. Among the 101 stations, 
15 (denoted by circles in the map on the right and by the dashed part of the right-hand 
histogram) are excluded from the analysis for validation. Dashed lines in the maps delimit 
the northern and southern slopes of the Alps. 
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and the Swiss Federal Office for Meteorology and Climatology (MeteoSwiss) . 
Annual maxima are extracted from daily snow depth measurements, which 
are read off a measuring stake at around 7.30 AM daily from November 1st 
to April 30th, for the 43 winters 1965-1966 to 2007-2008; we use the term 
"winter 1966" for the months November 1965 to April 1966, and so forth. 
Examples of such time series can be found in the Supplementary Materials, 
Blanchet and Davison (2011). As Figure 1 shows, the stations are denser 
in the Alpine part of the country, which has high tourist infrastructure and 
increased population density and traffic during the winter months. Their 
elevations range from 250 m to 2500 m above mean sea level, with only 
two stations above 2000 m. In order to validate our final model, we used 
86 stations to choose and fit the model and retained 15 stations for model 
validation. 

2.2. Marginal analysis and transformation. Let Z(x) denote the annual 
maximum snow depth at station x of the set X, which here denotes Switzer- 
land. Data are only available at the stations x £ D C X , so modeling Z(x) 
involves inference for the joint distribution of {Z{x),x G X} based on obser- 
vations from D, and extrapolation to the whole of X. In particular, as the 
station elevations lie mainly below 2000 m, any results must be extrapolated 
to elevations higher than 2000 m. 

Daily snow depths at a given location x are obviously temporally depen- 
dent. However, time series analysis suggests that, for every location x (zT) 
and every winter, daily snow depths show only short-range dependence. 
Hence, distant maxima of daily snow depths seem to be near-independent 
and, therefore, the D{un) condition for independence of extremes that are 
well separated in time [Leadbetter et al. (1983), Section 3.2] should be sat- 
isfied. Extreme value theory is then expected to apply to annual maximum 
snow depth: Z{x) at a location x may be expected to follow a generalized 
extreme- value (GEV) distribution [Coles (2001)] 



where = max(n, 0) and /i(x), a{x) > and ^(x) are, respectively, location, 
scale and shape parameters. 

Characterizing the probability distribution of Z{x) for all x X is equiv- 
alent to characterizing the probability distribution of f{Z{x)} for any bi- 
jective function /, which may be easier for a well-chosen /. A first step in 
our analysis is to transform the data at the stations to the unit Prechet 
scale. Whatever the values of the GEV parameters /x(x), a{x) and (,{x), 
taking f{z) = — l/logG(2) transforms {Z{x),x G X} into a spatial pro- 
cess {Z*{x),x G X} having unit Frechet marginal distributions, G*{z) = 
exp(— l/z). As it is easier to deal with Z* in general discussion, we will as- 
sume below that the time series at each station has been transformed in this 
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way. To do so, one might model the GEV parameters n{x), a{x) and ^{x) 
as smooth functions of covariates indexed by x, such as longitude, latitude 
and elevation [Padoan, Ribatet and Sisson (2010)]. However, due to the 
very rough topography of Switzerland and the influence of meteorological 
variables such as wind and temperature, snow depth exhibits strong local 
variation and additional covariates are necessary. A systematic discussion of 
such covariates and associated smoothing is given by Blanchet and Lehning 
(2010). The focus in the present paper is spatial dependence, so rather than 
adopt their approach, here we simply use GEV fits for the individual sta- 
tions to transform Z{x) at station x ^T) into Z*[x). Diagnostic tools such 
as QQ-plots showed a good fit even at low altitudes. 

2.3. Spatial dependence and regional patterns. A simple measure of the 
dependence of spatial maxima at two stations x,x' ^ X is the extremal co- 
efficient Oxx'- If Z*[x) is the limiting process of maxima with unit Frechet 
margins, then [Coles (2001), Chapter 5] 

(2) pr{Z*(x)<z,Z*(2;')<^} = exp(-6l^^//z), z > 0. 

One interpretation of Oxx' appears on noting that 

pr{Z*(a;') > z\Z*{x) > z]^2- O^x', z^oo. 

If ^xx' = 1; then the maxima at the two locations are perfectly dependent, 
whereas if Oxx' = 2, they are asymptotically independent as z ^ oo, so very 
rare events appear independently at the two locations. Although they do not 
fully characterize dependence, such coefficients are useful summaries of the 
multidimensional extremal distribution. In particular, it may be informative 
to compute all extremal coefficients {9xx'-,x' G X} for a given station x to 
see how extremal dependence varies. Figure 2 depicts such maps for the 
snow depth data, for four different reference stations x. Extremal coefficients 
{Oxx'^x' € 15} were estimated by the madogram-based estimator of Cooley, 
Naveau and Poncet (2006), and then kriged to the entire area using a linear 
trend on absolute altitude difference between x and x' . Similar maps have 
been proposed for gridded data by Coelho et al. (2008). 

Much information can be gleaned from Figure 2. A strong elevation effect 
is clearly visible. The map for Adelboden also suggests a directional effect: 
for this mid-altitude station in the Alps, there is more dependence with other 
middle-altitude stations in a roughly north-easterly direction. Another strik- 
ing feature visible in the two lower maps is near-independence between the 
northern and southern slopes of the Alps. Further such maps suggest the 
presence of the two weakly dependent regions separated by the black dotted 
line in Figure 1. A similar north/south separation was seen in Blanchet, 
Marty and Lehning (2009), for good reason: extreme snowfall events occur- 
ring in these two regions typically do not stem from the same precipitation 
systems. Whereas extreme snowfall events on the northern slope of the Alps 
usually arise from northerly or westerly airflows [Schiiepp (1978)], those in 
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Fig. 2. Extremal coefficient computed relative to Koppigen, Adelboden, Davos and Maloja 
(white points), estimated by the Cooley, Naveau and Poncet (2006) madogram estimator, 
and then kriged to the whole of Switzerland using a linear trend on absolute altitude dif- 
ference. 

the southern slope usually come from the south or south-west. These are less 
frequent, but when they occur they can be very severe, due to the proximity 
of the Mediterranean Sea. As snow cover results from the accumulation of 
many snowfall events during the winter, one can expect annual maximum 
snow depths on the northern and southern slopes of the Alps to be some- 
what disconnected. The winter of 1981 illustrates this: little snow fell on the 
southern slope of the Alps, while the northern slope received large amounts. 
Figure 2 nevertheless suggests that these two regions are asymptotically 
weakly dependent, since O^x' is generally larger than 1.7, but not necessarily 
asymptotically independent. Even between well-separated stations, Oxx' is 
rarely very close to 2, perhaps owing to the rather small area under study, 
in which the largest distance between stations is around 350 km. 

3. Spatial maxima. 

3.1. Max-stable processes. The spatial dependence highlighted in Sec- 
tion 2.3 suggests that we model Z*{x) as a spatial process of extremes. 
A max-stable process with unit Frechet margins is a stochastic process 
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{Z*{x),x € X} with the property that, if Z^*-^^(x), . . . , Z^*^^(x) are n inde- 
pendent copies of the process, then [de Haan (1984)] 

< max Zf-\{x),x G > has the same distribution as {nZ*{x),x G X}. 

li=l,...,n ^ ' J 

A consequence of this definition is that all finite-dimensional marginal dis- 
tributions are max-stable: if {xi, . . . ,xd} is a finite subset of X, then for all 

n G N, 

pr{Z*(a;i) < nzi, Z*{xd) < nzoY 

= pr{Z*(a;i) < zi, . ..,Z*{xd) < zd}, zi,...,zd>0. 
Such processes have several representations, two of which we now sketch. 

3.2. Smith's storm model. A general method of constructing max-stable 
processes is due to de Haan (1984). Let {{r]i,Si),i G N} denote the points of 
a Poisson process on (0, oo) x S with intensity r]~'^dr] x v{ds), where S is an 
arbitrary measurable set and v \s a. positive measure on S. Let {/(s,x),s G 
5,x G denote a nonnegative function for which, for all x ^X, 

f{s,x)u{ds) = 1. 

es 

Then the random process 

(3) Z* = \mayi{r]if{si,x)},xeX\ 

is max-stable with unit Prechet margins. Smith (1990) gives a rainfall-storms 
interpretation of this construction. He suggests regarding 5 as a space of 
storm centers, of /(s, •) as the shape of a storm centered at s, and of r] as 
a storm magnitude. Then r}f(s,x) represents the amount of rainfall received 
at location x for a storm of magnitude rj centered at s and Z*(x) in (3) is 
the maximum rainfall received at x over an infinite number of independent 
storms. 

Additional assumptions are needed to get useful models from (3). Smith 
(1990) proposes taking S = X = R^, letting u be the Lebesgue measure and 
f{s,-) be a multivariate normal density with mean s and covariance ma- 
trix S, that is, 

/(S, X) = (2^)~-°/2|5.|-l/2 gxp|_ 1 _ 5)^5^-1(2. _ g)}^ g g 

The resulting bivariate distribution of Z* defined by (3) at two stations xi 
and X2 is then 

pv{Z*{xi)<Zl,Z*ix2)<Z2} 

(4) 

f ^ ^fo. 1, Z2\ 1, zi 

= exp\ $ - + -log— ^[ - + -log — 

I zi \ 2 a z\) Z2 \1 a Z2 



8 J. BLANCHET AND A. C. DAVISON 



Isotropic case Anisotropic case 




Fig. 3. Smith's process in two dimensions with two different matrices E — (Ttjtj/)tj_^/g{i_2} • 
Upper left image: a simulated field with th = T22 = 17 and ri2 = (isotropic case). Upper 
right image: a simulated field with th — 25^, T22 = IS'^ and T12 = 14'^ (anisotropic case). 
Lower images: corresponding pairwise extremal coefficient. 



where $ is the standard normal distribution function and a is the Mahalano- 
bis distance given by 

(5) = (xi - X2)^S~"^(X1 - X2). 

Below we will call this model the Smith process. 

Two simulated Smith processes with different matrices S are shown in 
the top row of Figure 3. The anisotropic case arises when S is not spherical, 
that is, not of the form S = t'^Id, where > and Id is the identity matrix 
of side D. The resulting geometric anisotropy [e.g., Journel and Huijbregts 
(1978)] can easily be seen by computing pairwise extremal coefficients. Tak- 
ing zi = Z2 = z in (4) gives, according to (2), 

(6) 9,,,, = 2<P{a/2). 

The Mahalanobis distance a appearing in (6) gives different weights to the 
different components of the vector (xi — X2). The limiting cases a ^ O"*" and 
a— >■ +00 correspond, respectively, to perfect dependence, 6x1x2 = 1) a-iid hi- 
dependence, 6x1x2 = 2. For a given station xi, surfaces {x2 G X, 6x^^x2 = c} 
are, according to (6), such that (5) is constant. If T, is spherical, then such 
surfaces are circles in two dimensions and spheres in three dimensions. Oth- 
erwise, they are ellipses and ellipsoids, respectively. 
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3.3. Schlather's storm model. A second method of construction of max- 
stable processes was proposed by Schlather (2002). Let {r/j, z G N} denote the 
points of a Poisson process on with intensity •q~'^drj. Let {W{x),x € X} 
be a stationary nonnegative process that satisfies E[VF(x)] = 1 for all x € X, 
and let Wi, z € N, be independent copies of this process. Then [Schlather 
(2002)] the random process 

(7) Z* = {ma^r]iWi(x),x£x} 

is max-stable with unit Prechet margins. When Wi{x) = f{x — Si), where / is 
a density function on X and the Si are the points of a Poisson process with 
unit rate on a measurable set S, then (7) is equivalent to the storm model of 
Section 3.2. Smith's model (4) corresponds to taking / to be a multivariate 
normal density, extended by de Haan and Pereira (2006) to Student t and 
Laplace densities. Like Smith's model, the model (7) has a simple interpre- 
tation: the rjW are spatial events all having the same stochastic dependence 
structure but differing in their magnitudes rj. An appealing difference be- 
tween this and the Smith model is that the shapes of the events may vary 
if the process W permits this. 

Additional assumptions are again needed to get useful models from (7). 
Schlather (2002) proposes taking W to be the positive part of a station- 
ary Gaussian process with correlation function p, scaled so that E[max{0, 
Ty(x)}] = 1 for all x X. He shows that the corresponding bivariate distri- 
bution of Z* at two stations xi and X2 is 

pr{Z*(xi) < zi,Z*{x2) < Z2} 

(8) 

= exp|-V- + -Vl + ./l-2(p(/.) + l): 



2v^i ^ayV V (2:1+22)^, 
where h € M_|_ is the Euclidean distance \\x2 — xi|| between the two stations. 
Below we call this max-stable model Schlather's process. 

A simulation from an isotropic version of this model with X corresponding 
to Switzerland is shown in Figure 4. The isotropy can be easily seen by 
computing pairwise extremal coefficients. Taking zi = Z2 = z in (8) gives, 
according to (2), 

1/2 



(9) 0X1X2 = 1 + 



1 - p(\\xi - X2 



Here the extremal coefficients involve the Euclidean distance between the 
two locations. For a given station xi, surfaces with the same extremal coef- 
ficents c G [1,2], that is, surfaces {x2 € X, 6x^x2 = c}, are, according to (9), 
such that ||xi — X2II = c' . Such surfaces are circles in two dimensions and 
spheres in three dimensions. The limiting case ||xi — X2II ^ O"*" corresponds 
to perfect dependence, 9x^^x2 = 1- If, like most geostatistical correlation func- 
tions, the underlying Gaussian process has p{h) when h — > +00, then the 
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Fig. 4. Schlather's process in two dimensions for a Cauchy covariance function 
p{h) = (1 + /19^)~^ . Left image: one simulated field. Right image: corresponding pair- 
wise extremal coefficient. 

limiting case \\xi — X2\\ +oo corresponds to 9xiX2 = 1 + 2~^/^ ~ 1.707, and 
so independent extremes do not arise even at very large distances. Moreover, 
as an isotropic correlation function can give correlations no smaller than 
-0.403 in M2 and -0.218 in [Matern (1986), page 16], under Schlather's 
model we have 9x^x2 ^ 1-838 for any xi, X2 in and 6x1x2 ^ 1-780 for any 
xi, X2 in R^. Thus, it is impossible to produce independent extremes using 
such a process, no matter how distant the stations. Davison and Gholam- 
rezaee (2010) have proposed extensions to allow independence in (8), and 
Kabluchko, Schlather and de Haan (2009) have extended both Smith's and 
Schlather's representations. 

4. Max-stable process for extreme snow depth. 

4.1. General. As pointed out in Section 2.3, snow depth data show two 
key characteristics that should be explicitly modeled in the max-stable pro- 
cess. First, dependence is anisotropic, due to the strong elevation effect and 
the presence of a main direction of dependence. Second, Switzerland seems 
to be divided into two weakly dependent climatic regions: the northern slope 
of the Alps together with the Plateau, which is the low altitude region north 
of the Swiss Alps; and the southern slope of the Alps. In this section we 
propose to extend the Smith and Schlather models of Sections 3.2-3.3 to 
account for these features. Other representations described in Kabluchko, 
Schlather and de Haan (2009), in Davison and Gholamrezaee (2010) or in 
Davison, Padoan and Ribatet (2010) are not considered in this paper. 

4.2. Modeling anisotropy. Smith's model can directly model anisotropy 
using a nonspherical T, matrix in (4). The simple version of Schlather's 
model is isotropic, but it can easily account for anisotropy by considering 
a transformed space <f instead of X. 
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Anisotropy of Smith's model arises from the fact that the distance used 
in the extremal coefficient (6) is Mahalanobis distance (5) rather than Eu- 
clidean distance. Using the eigendecomposition S = UAU'^ , where U is a ro- 
tation matrix and A a diagonal matrix of positive eigenvalues, we may write 



roots of the diagonal elements of A. If Ai denotes the first element of A, 
then (10) can be written as = X^W^V, where V = Xl^"^ A^^^'^U . The 
squared Mahalanobis distance (5) is 



which is exactly that between xi = Vxi and X2 = Vx2 in the isotropic case, 
that is, when using a D-dimensional spherical covariance matrix Xilo in (4). 
Thus, the anisotropic Smith model on X is just the isotropic Smith model 
on the transformed space 

Similar ideas can be used with Schlather's model, by applying it on ^ = 
VX, where in three dimensions we may take 



as for Smith's model. In the rest of the paper we will use the term climate 
space for the transformed space $ = VX in which isotropy is achieved. Fig- 
ure 5 illustrates the climate space transformation, allowing an anisotropic 
Schlather model, with the same V matrix as that corresponding to the 
anisotropic case of Figure 3. Compared to Figure 4, constant extremal coef- 
ficients correspond to ellipses, allowing us to model directional effects. 





(11) 




■+' 




Fig. 5. Anisotropic Schlather model resulting from climate space transformation. Left 
image: a simulated field. Right image: corresponding extremal coefficients. 
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Geometric anisotropy as induced by the V matrix is a special case of ran- 
ge anisotropy [Zimmerman (1993)]. In the nonextremal framework, this 
idea has been extended to nongeometric range anisotropic models, in which 
nested covariances are used with different range parameters in different di- 
rections, but in general this does not define a valid covariance function. 
Ecker and Gelfand (2003) introduced product geometric anisotropy, under 
which covariance functions are products of geometric anisotropic covari- 
ances. Space transformation has also been used by Sampson and Guttorp 
(1992) to model nonstationary spatial covariance structures, allowing more 
complex transformations than the affine transformation considered here. In 
addition to these global methods, local methods for modeling anisotropy 
and more general forms of nonstationarity also exist. These can be divided 
in three main families [Schabenberger and Gotway (2005)]. The moving win- 
dow approach of Haas (1990) estimates a covariance function locally within 
a neighborhood. The convolution method of Higdon (1998) allows the con- 
struction of weakly nonstationary processes by convolving a zero-mean white 
noise process with a kernel function whose parameters can depend on loca- 
tion. The method of weighted stationary processes [Fuentes (2001)] allows 
one to write the nonstationary covariance function as a weighted mixture 
of isotropic covariances, where the weights depend on the location. Fuentes, 
Henry and Reich (2010) use a Dirichlet process mixture as the basis for 
a flexible copula approach to space-time modeling of extreme temperatures, 
but it does not correspond to a max-stable process model, and the relatively 
long-range dependence of temperatures can be modeled more simply than 
can precipitation phenomena such as rain- and snowfall. It would be very 
valuable to apply these ideas in the max-stable context, but the unavailabil- 
ity of a likelihood function seems to be a major obstacle. 

The idea of space transformation was used by Cooley, Nychka and Naveau 
(2007) in modeling US precipitation. Instead of using the three-dimensional 
geographical coordinates (longitude, latitude, elevation) for locating sta- 
tions, the authors work in a "climate space," namely, the two-dimensional 
space given by elevation and mean precipitation for the months April to 
October. Unlike in Cooley, Nychka and Naveau (2007), our transformation 
is affine, giving more weight to elevation through C3, and defining a main 
direction of dependence along the o-axis. A higher-dimensional space could 
of course be used for X. In particular, one could use the four-dimensional 
space of (longitude, latitude, elevation, mean snow depth), thus blending 
the Cooley, Nychka and Naveau (2007) approach with ours; see Section 6. 

4.3. Modeling climate regions. Different approaches to accounting for 
the impact of the climate regions on the extremes are possible: 

1. The climate regions are independent. This is equivalent to saying that 
two max-stable processes govern the two regions independently. In terms 
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of spatial dependence, extremal coefficient maps will be of the form of 
Figure 3 or 5 but replacing the Swiss border by the border of the northern 
region alone for the pairwise dependence with a station located in the 
north, and similarly for the southern region. 

The climate regions are weakly dependent. Since dependence between pair 
of stations decreases when distance increases, one way to model weakly 
dependent regions is to increase the distance between them. This can be 
done by adding to X a coordinate equal to in the northern region, and to 
1 in the southern region. If the other coordinates are (longitude, latitude, 
elevation), then the V matrix of the climate space transformation (11) can 
be written in the most general case as a 4 x 4 matrix with one column 
comprising apart from one element. Nevertheless, for computational 
reasons it may be better to consider the rotation matrix U of Section 4.2 
as being a rotation matrix in the (longitude, latitude) plane and thus to 
set 

/ cos a — sin a \ 
-J , _ C2sina C2Cosa 1 
C3 ' 

\ C4/ 

we shall do this in Section 6. In the four-dimensional climate space 
^ = VX, the squared distance between two stations xi and X2 is {l^(a;i — 
^2)}'^ {y {"^1 — X2)}- But the fourth coordinate of xi — 3:2 is if the two 
stations belong to the same region and ±1 otherwise. The squared dis- 
tance will then equal that in the (longitude, latitude, elevation) climate 
space if the two stations are in the same region, and be increased by c| 
otherwise. We thus increase the distance between the climate regions, and 
therefore decrease the dependence between them, without increasing the 
distance between stations of the same region. To see how the extremal 
coefficients behave, see the left-hand side of Figure 6. 
The climate regions are weakly dependent in continuous space. Since the 
additional coordinate introduced above jumps from to 1 at the border 
between the regions, it induces a discontinuity of the extremal coefficients 
which is visible in the left map of Figure 6; see the cyan and magenta 
ellipses. This seems unrealistic and something smoother is preferable. An 
easy way to impose space continuity is to take the border to be a band in- 
side which the fourth coordinate is linearly interpolated between and 1, 
with value on the upper-border of the band and 1 on the lower-border. 
With this simple interpolation, there is no jump at the border and curves 
of constant extremal coefficient are continuous, as in the right-hand side 
of Figure 6. The width of the band must be estimated from the data; we 
return to this in Section 5. 
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Discontinuous space Continuous space 




Fig. 6. Example of extremal coefficients with weakly dependent regions in discontinuous 
space (left image) and continuous space (right image). The images are the same, except 
in a 10 km wide hand around the north/south border (dashed line). 

5. Model estimation and selection. 

5.1. Pairwise likelihood. Statistical inference for parametric models is 
ideally performed using the likelihood function. Let D= {xi, . . . ,X£)} C 
X denote the 86 stations whose maxima are used for fitting the mod- 
els. Computation of the likelihood requires the joint density function of 
{Z*{xi), . . . ,Z*{xd)}, but in the framework of max-stable processes, this 
is infeasible because only the bivariate marginal distributions are available. 
Padoan, Ribatet and Sisson (2010) proposed replacing the full likelihood by 
a pairwise likelihood function [Cox and Reid (2004); Varin (2008)]. This idea 
is also used by Davison and Gholamrezaee (2010), Davison, Padoan and Ri- 
batet (2010) and by Smith and Stephenson (2009), the latter in a Bayesian 
framework. 

Let Zik denote the kth. observed maximum for the ith station, trans- 
formed so that time-series (zji, . . . , Zjx) at each station have unit Frechet 
distributions; here k £ {1, . . . , K}, with K = 43 years, and i € {1, . . . , D}, 
with D = 86 stations. Let (3 = {(3i, . . . ,/3r) denote the parameters to be es- 
timated. Then the pairwise marginal log-likelihood is 

K 

(13) = Y,Y.^og f{zik,z,k;P), 

k=l i<j 

where /(•,•) is the bivariate density of the unit Frechet max-stable pro- 
cess, that is, the derivative of equation (4) for Smith's model or of (8) for 
Schlather's model, and the second summation is over all distinct pairs of 
stations, D(^D — l)/2 terms in all. Under suitable regularity conditions, the 
maximum pairwise maximum likelihood estimator /3 has a limiting normal 
distribution as K +oo, with mean /3 and covariance matrix of sandwich 
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form estimable by H{/3) ^ J{j3)H{$) ^, where 

,14, = 2 

k=l i<j ^ ^ 

j(a\ _ ^ s:^ 9 log f{zik,Zjk;f3) dlogf{zik,Zjk;f3) 

(15) APJ-l^l^ 

k=l i<j ^ ^ 

are the observed information matrix and the squared score statistic cor- 
responding to £p. The use of the pairwise hkehhood estimator for Smith's 
process was validated by Padoan, Ribatet and Sisson (2010) in a simulation 
study. 

5.2. Estimation in practice. Estimating the maximum pairwise likeli- 
hood estimator requires the maximization of (13) with respect to the R pa- 
rameters. We found that the R function optim gave quite poor results for 
our application: the surface can have many local maxima, and optim 
and similar functions find it hard to deal with them. After some experi- 
mentation, we therefore adopted a profile likelihood method. Given a set of 
(i? — 1) parameters it is easy to find the value of j3r that maximizes 
the single- variable function ^p(-,/3_f.). This suggests the following iterative 
algorithm: 

1. Take initial parameters (3 = (/3i, . . . , /3r). 

2. For r in 1, . . . , i2: 

(a) find the value (3r that maximizes the pairwise likelihood with respect 
to the scalar (3^-, holding the other parameters, fixed, that is, 

/3r = arg max £p (/3r , /3-r ) ; 

(b) then update the rth component of (3 to $r- 

3. Go to step 2, stopping when no change to any Pr can increase the pairwise 
log- likelihood. 

To assess the performance of this algorithm, we simulated 200 data sets, 
each comprising 43 independent copies of Schlather's max-stable random 
field (8) with Cauchy covariance function p in a three-dimensional climate 
space. Each of the copies is observed at the same D = 100 stations, so the 
number of observations is very similar to those for the annual snow depth 
data; see Section 2.1. The climate transformation is defined through a 3 x 3 
matrix V as in (11). The model has three parameters for the V matrix and 
two for the covariance function, which induces middling dependence: about 
25% of the pairs of stations have extremal coefficients 9 < 1.68; recall from 
Section 3.3 that for Schlather's model, 9 < 1.707. We started from the same 
initial point for each of the 200 data sets and optimized the log pairwise 
likelihood (13) with (i) eight optimization procedures within the R function 
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Fig. 7. Difference in pairwise log-Ukelihood at convergence between the profiling algo- 
rithm and eight algorithms for simultaneous parameter estimation, for 200 simulated data 
sets. The eight algorithms are: Nelder-Mead, NM; the quasi- Newton method of Broyden, 
Fletcher, Goldfarb and Shanno, BFGS; three conjugate gradient methods, CG-FR, CG-PR 
and CG-BS; and a variant of simulated annealing using starting temperatures of 5, 10 and 
20, SANN-5, SANN-10 and SANN-20. The help for the R function optim gives more details 
of these algorithms. 



optim with all parameters estimated jointly [Blanchet and Davison (2011)]; 
and (ii) the above profile likelihood algorithm. Figure 7 shows the differences 
between the pairwise log-likelihoods for the methods at convergence, for the 
200 data sets. The profiling method never gives lower maximized pairwise 
likelihoods than the other algorithms, and they are almost always higher. 
Further simulations with small- and large-range dependence gave similar 
results [Blanchet and Davison (2011)]: overall profiling is clearly better than 
the other algorithms. Those that compare best with profiling, viz., Nelder- 
Mead and simulated annealing, are designed for rather rough surfaces with 
many local optima. These simulated data are relatively simple compared to 
the real data, which are neither exactly unit Frechet after transformation 
from (1) nor follow a pure max-stable process. Furthermore, the max-stable 
model used for the simulation is quite simple, with only five parameters to be 
estimated, so the profiling approach seems necessary for our, more complex, 
application. 

5.3. Model selection. Model selection criteria play an important role in 
deciding which of the fitted models should be preferred. As in Padoan, Rib- 
atet and Sisson (2010), we propose to use the composite likelihood informa- 
tion criterion [Varin and Vidoni (2005)], which extends the TIC [Takeuchi 
(1976)] to the composite likelihood setting, and is defined as 

CLIC = -2ip{$) + 2tr{H{$)-^J0)}, 

where H and J are, respectively, the observed information matrix and 
the squared score statistic corresponding to ip, defined at equations (14) 
and (15), and $ is the maximum pairwise likelihood estimator. Lower values 
of CLIC correspond to better quality models. 



SPATIAL MODELING OF EXTREME SNOW DEPTH 



17 



6. Application to snow depth in Switzerland. 

6.1. Fitted models. We fitted the different models described in Section 4 
to our snow depth data, using both Smith and Schlather max-stable struc- 
tures for the extremes. For Schlather's model, different choices of Gaus- 
sian covariance function p lead to different distributions (8). We used nine 
such functions, namely, the spherical, circular, cubic, Gneiting, exponential, 
Matern, Gaussian, powered-exponential and Cauchy covariance functions 
[Banerjee, Carlin and Gelfand (2003); Schabenberger and Gotway (2005)]. 
Each has either one or two parameters and the first four have an upper 
bound. They all are such that p{h) 1 when /i — )• 0"^ and p{h) — )• when 
h +00. As mentioned in Section 3.3, this constrains the extremal coeffi- 
cient for Schlather's model to correspond to dependent data. Nevertheless, 
we will see that such an assumption seems justified in our case. 

The coordinates x we considered are geographical coordinates (longitude, 
latitude, elevation), region number (see Section 4.3) and mean snow depth 
during the winters 1966-2008. Mean precipitation was considered as a pos- 
sible climate coordinate in Cooley, Nychka and Naveau (2007) 's study of 
extreme precipitation. The idea of using mean snow depth is that stations 
with similar snow depth are probably influenced by the same weather pat- 
terns and should therefore be closer in the climate space than are stations 
with different snow cover. Other climate variables that could be considered 
are temperature, wind direction and wind speed, which are also measured 
at the stations, but these values are of relatively poor quality with many 
missing values, so we decided not to use them. 

In addition to the models illustrated in Section 4, we allowed the possi- 
bility of having different climate spaces in northern and southern regions, 
that is, to have different climate space transformation matrices V. In three 
dimensions, for example, two V matrices as in (11) will have to be esti- 
mated, with a total of 6 parameters. In the continuous-space case illustrated 
in Figure 6, all coefficients a and c are linearly interpolated around the 
north/south border. We also considered different mixtures of the above pos- 
sible coordinates. In all cases, we used longitude and latitude, plus possibly 
the elevation, region number and mean snow depth, or combinations of these 
three coordinates. In total, 65 types of models were considered, each of them 
being estimated for one Smith and nine Schlather processes, giving 650 fits 
in all. A description of the 65 model types is given in the Supplementary Ma- 
terials [Blanchet and Davison (2011)]. All were estimated using the iterative 
profiling algorithm of Section 5.2. 

6.2. Model comparison. A summary of the CLIC values for the 585 fit- 
ted Schlather models, rescaled by division by D — 1 in order to give log- 
likelihood values that would correspond to independent data, is shown in 
Figure 8. There are relatively small differences among them, though the 
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Fig. 8. Rescaled CLIC values for all 65 x 9 fitted Schlather models. 

Gneiting and Gaussian covariance functions seem to perform less well and 
the spherical and circular covariance functions have the 25 best CLIC val- 
ues. These covariance functions have an upper bound and are governed by 
only one parameter. Schlather's model always performs better than Smith's 
model, whatever the chosen covariance function: the rescaled CLIC with 
Smith's model is between 30 and 300 units higher than with Schlather's 
model, with a minimum value of 15,650 attained for model 47. Whether 
with Smith or Schlather models, the same patterns appear. In particular, 
the first eight models, which perform poorly, correspond to models in Eu- 
clidean space, without climate space transformation. The benefit of working 
in a transformed space in order to allow for anisotropy is thus clear. This 
effect is particularly striking for Smith's model, for which it is equivalent 
to saying that a nonspherical S matrix (see Section 3.2) should be used: 
there is a difference of 300 between the lowest rescaled CLIC values in the 
Euclidean and climate spaces. The models numbered 10, 11, 17, 21, 25, 
29, 30, 36, 40, 44, 45, 51, 55, 59 and 63, which are also poor, correspond 
to cases when neither elevation nor the mean snow depth are considered 
[Blanchet and Davison (2011)]. As the mean snow depth is strongly related 
to elevation, the latter is a very important climate coordinate. It seems to 
be more informative than the mean snow depth; models using elevation but 
not mean snow depth as a coordinate always have lower CLIC values than 
in the converse case. 

6.3. Selected model. According to the CLIC, the best fit is given by 
Schlather's model with spherical covariance function, and a 5-dimensional 
climate space X of coordinates (longitude, latitude, elevation, region num- 
ber, mean snow depth) with different transformations in the north and south 
but imposing space continuity; this, model number 47 in Blanchet and Davi- 
son (2011) has a CLIC = 15,611.66. This means that two V matrices are 
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Table 1 

Parameters (standard errors) of the model selected by CLIC; Schlather's model with 
spherical covariance function, two climate transformations but a continuous space (the 
band around the north/south border is about 5 km wide) 

Covariance parameter 

447.45 (43.32) 



Climate space parameters 



Main direction Latitude (km) Elevation (km) Mean snow Region num. 

(radian) depth (cm) 

North 0.36 (0.03) 4.98 (0.80) 274.7 (35.4) 1.26 (0.46) x 

South 0.17 (0.06) 4.70 (1.16) 406.5 (161.3) 6.41 (2.45) 449.4 (37.4) 



estimated, each of the form (12) but in five dimensions, and thus having five 
parameters: the main direction of dependence, and the four parameters c 
associated to the latitude, elevation, region number and mean snow depth. 
Since the region number is a binary variable, the c value for the northern 
region can be fixed equal to zero. The range parameter of the spherical co- 
variance function and the width of the band between the regions are also 
estimated, for a total of 11 parameters, whose estimates and standard errors 
are shown in Table 1. As the pairwise likelihood is not differentiable with 
respect to the band width, no standard error is given for it. The second- and 
third-best fits are also obtained with spherical covariance functions with sim- 
ilar models as in Table 1 but without the mean coordinate (model number 49, 
with CLIC = 15,612.03) or the region number coordinate (model number 46, 
with CLIC = 15,612.56), that is, using a four-dimensional space X. In the 
latter case, values of the estimated coefficients are such that the northern 
and southern regions are disjoint in the climate space, although no region 
number coordinate is used to separate them. These two models perform 
similarly because the mean coordinate should provide information about 
the local variability of snow depth, part of which agrees with the regional 
division between the northern and southern slopes; thus, the mean coordi- 
nate and region number carry similar information. According to Figure 8, 
it seems better to use both coordinates, but using one of them increases the 
CLIC only very slightly. 

It is no surprise that in Table 1, elevation is the most influential coordi- 
nate in the climate distance, and thus in the dependence function. In the 
north, for example, dependence between two stations at the same eleva- 
tion but 10 km apart along the main direction of dependence, an angle of 
a = 0.36 radians in the sense of an Argand diagram, at the same elevation 
but 2 km apart perpendicularly to the main direction of dependence, and at 
the same latitude and longitude but 40 m apart in elevation, are all equal. 
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An interesting feature is the main direction of dependence in the northern 
region, which can be explained by two facts: 

1. due to the strong elevation effect, the north slope of the Alps (the moun- 
tainous part of the northern region) is very weakly dependent on the 
Plateau (the low-elevation part of the northern region). But both subre- 
gions are oriented along the North Alpine ridge, and dependence is thus 
higher in this direction; 

2. this direction is also broadly that of the two widest valleys in Switzerland, 
the Rhone and Rhine valleys, as shown by the main green valleys in Figu- 
re 1. These are wide enough to direct snow-bearing clouds along them, 
thus inducing strong directional dependence of precipitation. 

The high value associated to the region number coordinate gives the lowest 
possible dependence, 6xx' = 1-707, between extremes in the northern and 
southern regions. 

Figure 9 shows maps of the estimated pairwise dependence under the 
max-stable model of Table 1, obtained by extrapolating the mean snow 
depth at ungauged stations where no data are available. To do this, we 
performed spatial kriging with a spline dependence on elevation, to allow for 
the fact that temperatures at stations below 800 m may exceed 0°C even 
when it is snowing at higher altitudes, leading them to suffer rain rather 
than snow. The resulting smooth mean process was successfully validated 
on the additional 15 stations [Blanchet and Davison (2011)]. Figure 9 clearly 
shows both the elevation effect and the weak north/south dependence. The 
low bandwidth, of about 5 km, induces an abrupt change of the extremal 
coefficient around the north/south border. 

6.4. Model checking. For a first check on the quality of the selected 
model, we compare its predicted extremal coefficients, obtained by replac- 
ing the parameters involved in (9) by their estimates from Table 1, with the 
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Fig. 9. Pairwise extremal coefficient with Koppigen and Davos (white circles) predicted 
by the selected max-stable model. 
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Fig. 10. Extremal coefficient for pairs of stations as a function of the distance between 
them, in Euclidean space (left plot) or climate space (center and right). The red curve is 
the extremal coefficient curve for the corresponding max-stable model. 



naive estimators of Schlather and Tawn (2003) or the madogram-based esti- 
mator of Cooley, Naveau and Poncet (2006). As the extremal coefficients (9) 
are functions of distance between stations, we plot naive and predicted ex- 
tremal coefficients against distance. Figure 10 shows such comparisons for 
our selected model and for the best Smith model. For clarity, we only show 
the madogram-based estimator of Cooley, Naveau and Poncet (2006), with 
and without binning. The naive estimators of Schlather and Tawn (2003) 
give essentially the same picture, but with slightly higher variability. 

Figure 10 shows that the Smith model fits the data less well than the 
Schlather model. In particular, the extremal coefficient curve of the Smith 
model crosses the point cloud for the binned madogram, whereas our selected 
model follows it quite well up to a climate distance of 400, and then under- 
estimates it. A limit of about 1.8 would be expected from the madogram, 
but cannot be attained with Schlather's model; see Section 7. 

Another way to check our model is to compare the empirical distribution 
of maxima of subsets of stations, that is, = max{Z* {xi),Xi G A}, with 
maxima predicted by the selected model. The distribution of under the 
selected model is known analytically only when A comprises two stations, 
but samples of Z^ can be simulated for any A. Since realizations z*^ of Z*^ 
are available for K = 43 years, one can compare the empirical quantiles of 
Z^ with the simulated ones. More precisely, given a subset A, we simula- 
te M independent series z^^™^ of length K, and thus obtain M replicates of 
the observed Frechet series z^. Ordered values of observed z^ can then be 

compared with ordered values of the graphical test of fit. Pointwise 

and overall confidence bands can also be derived from these simulations 
[Davison and Hinkley (1997), Section 4.2.4]. 

Figure 1 1 uses this approach to compare fitted and empirical distributions 
for different groups of three or four stations taken from the 15 not used to fit 
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Fig. 11. Comparison of empirical and model quantiles for annual maxima of groups of 
stations not used in the fitting. The stations used for each panel are shown in its map, and 
the envelopes are 95% pomtwise and overall confidence hands obtained from M = 5,000 
simulations. 



the model, some groups being tightly clustered, and others being dispersed. 
The fit seems to be broadly satisfactory in all cases. Even the dependence 
between stations whose climate distance is larger than 500 units seems to 
be well-modeled, despite the mismatch between the fitted and empirical 
pairwise extremal coefficients at such distances seen in Figure 10. 

6.5. Risk analysis. For risk management it is important to be able to 
assess how extreme events are likely to occur in the same year in differ- 
ent places. A first answer to this question can be obtained by computing 
probabilities of the form pr[{Z*(x) > z,x ^ A}] for a group of stations A 
and different high levels z. Figure 12 plots such probabilities for different 
groups A when z is the r-year return level of the unit Frechet distribution. 
By back-transformation from equation (1), this is equivalent to computing 
the joint survival distributions pr[{Z(a;) > RLr(x),x € A}\ where RLr(x) 
denotes the r-year return level at station x, that is, the probability that all 
stations in A receive more snow a given year than their r-year return level. 
Under independence, this probability equals r"'-^' for any possible set A^ 
where |^| is the number of stations in A, whereas it equals r~^ under full 
dependence. Figure 12 shows very good agreement between the observed and 
predicted distributions using the model, whereas the risk is underestimated 
under the hypothesis of independent stations and overestimated under the 
hypothesis of full dependence. The underestimation is more striking for quite 
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Fig. 12. Risk analysis of groupwise annual maxima: joint survival probability versus 
return period. In the right-hand of each panel the envelope is a 95% pointwise confidence 
band obtained from M — 5,000 simulations. Stations indicated in green were not used for 
fitting. 



dependent stations, such as those in the left-hand panel of Figure 12. When 
distance increases, the difference between the dependent and independent 
cases is less striking but our max-stable model fits better even for pairs of 
stations that are 980 climate distance units apart; this is almost the largest 
climate distance between pairs of stations. The right-hand panel corresponds 
to a group of seven stations in the eastern Plateau. Our model clearly gives 
more realistic risk probabilities than does the independence assumption. Ex- 
treme snow events in the low-elevation Plateau generally occur over a large 
region due to the easy weather circulation. A typical example is the ex- 
traordinary snowfall event that occurred on March 5th 2006 over the entire 
Plateau, with snow measurements of 54 cm at Zurich, 49 cm at Basel and 
60 cm at Sankt Gallen. This was the largest snow depth recorded since 1931 
[Zanini, Sutter and Gerstgrasser (2006)]. 

7. Discussion. The models discussed here are a step toward modeling 
spatial dependence of extreme snow depth. They are based on the Smith 
(1990) and Schlather (2002) max-stable representations, designed to model 
extreme snow depth explicitly. In particular, they can account in a flexible 
way for the presence of weakly dependent regions. They involve a climate 
transformation that enables the modeling of directional effects resulting from 
phenomena such as weather system movements. In the proposed methodol- 
ogy, model fitting is performed by using a profile-like method for maximizing 
the pairwise likelihood function, and model selection is performed using an 
information criterion. 

We applied this methodology to 86 stations with recorded snow depth 
maxima. Performance of the selected model at small and large scales was 
assessed on these stations, together with 15 other stations, by comparing 
empirical and predicted distributions of group of stations. By accounting 
for spatial dependence, our model gives clearly more realistic probabilities 
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of extreme co-occurrence than would a nonspatial model. Such quantities 
are important for adequate risk management. 

Considered as a whole, the max-stable models proposed in this paper con- 
stitute a family of flexible models that could potentially be applied to other 
kinds of climate data, in particular, extreme precipitation and temperature. 
Further improvements could nevertheless be investigated, as discussed be- 
low. 

In this paper we focus on modeling the spatial dependence of extremes, 
rather than on the marginal distributions. A first step was thus to transform 
maxima from their original scale to a common unit Frechet distribution. In 
the application to snow depth data, this transformation was done by us- 
ing the GEV distributions fitted to the time series, considered separately. 
A fuller spatial model would consider the three GEV marginal parameters 
as response surfaces. Using the models presented in this paper, one could 
then simultaneously estimate the spatial dependence and the spatial inten- 
sity of maxima, following Padoan, Ribatet and Sisson (2010) and Davison 
and Gholamrezaee (2010). These authors use simple functions of longitude, 
latitude and elevation, but the very complex Alpine topography results in 
an extremely variable pattern of snow, and we were unable to find satisfac- 
tory marginal response surfaces for our application. Blanchet and Lehning 
(2010) describe other approaches that appear to be more satisfactory, but 
modeling of the margins requires more investigation. Time could be used 
as a covariate in order to allow for the potential impact of climate change 
on extreme snow events; for example, the retreat of the glaciers is strongly 
affecting microclimates at high altitudes. This notwithstanding, exploratory 
work suggests that although climate change has affected mean snow lev- 
els [Marty (2008)], its effect on extreme snow events is not yet discernible, 
except possibly at low elevation [Laternser and Schneebeli (2003)]. 

A second improvement might be the consideration of event times, which 
could be incorporated into the pairwise maximum likelihood procedure [Ste- 
phenson and Tawn, (2005)]. For our data, the co-occurrence of annual max- 
ima is quite variable. For winters such as those of 1975 and 2006, snow depth 
reached its maximum almost simultaneously all over Switzerland. For win- 
ters such as those of 1980, 2007 and 2008, the annual maxima occurred at 
quite different dates; see the Supplementary Materials, Blanchet and Davi- 
son (2011). Including this information by modifying the pairwise likelihood 
contribution of maxima occurring simultaneously at two stations might yield 
more precise inferences, as shown in Davison and Gholamrezaee (2010). 

Last but not least, this article has used only snow data gathered from mea- 
surements in flat, open and not too exposed fields. Extrapolation to steep, 
windy and forest terrains may thus be unsatisfactory. In particular, prefer- 
ential deposition of snow [Lehning et al. (2008)] may imply that snow depth 
on slopes is more extreme than on representative flat fields. This could have 
important implications for avalanche risk [Lehning et al. (2006)] but could 
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not be considered here due to lack of data. This could be investigated us- 
ing data from automatic stations located at higher elevations, mostly above 
2,200 m, and in various terrains, though such data are unfortunately avail- 
able only for about ten years. A spatial model for exceedances over high 
thresholds [Davison and Smith (1990)] would be a valuable addition to the 
extreme- value toolkit for dealing with spatially-dependent short time series. 
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SUPPLEMENTARY MATERIAL 

Supplementary Material for "Spatial modeling of extreme snow depth" 

(DOI: 10.1214/11-AOAS464SUPP; .pdf). This contains example time series 
of data, and further discussion of the estimation algorithm and of the fitted 
models. 
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